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We present an analytic proof demonstrating the equivalence between the Random Phase Ap- 
proximation (RPA) to the ground state correlation energy and a ring-diagram simplification of the 
Coupled Cluster Doubles (CCD) equations. In the CCD framework, the RPA equations can be 
solved in 0(N 4 ) computational effort, where N is proportional to the number of basis functions. 

There has recently been a revival of interest in RPA in the quantum chemistry community. The RPA is popular 
for calculations of excitation energies both in finite systems*^ and in solids,—^ and is related to time-dependent 
density functional theory^&i As a technique for describing electronic correlations, RPA has significant advantages, 
particularly for those interested in density functional theory. It describes dispersion and van der Waals interactions 
correctly and is exact for long-range correlations^ Left-right static correlations seem to be properly described by 
RPApi and RPA fixes the pathologies of nonlocal Hartree-Fock-type exchange in metallic systems. Readers interested 
in details about RPA for ground state correlation can refer to the recent paper by Furcho^ where he discusses 
an interesting simplification to reduce the computational cost of RPA correlation and provides ample background 
information about RPA. Note that his work focuses on direct RPA, in which the exchange terms are neglected; as 
discussed later in this communication, this is the form of RPA most useful in the context of density functional theory. 

A connection between the RPA correlation energy and a ring diagram approximation to CCD was first mentioned by 
Freeman in his 1977 paper Very recently, A. Griineis and G. Kresse reproduced this evidence and found numerical 
proof of the equivalence between these two approaches.— Here, we offer an analytic proof that these two problems 
yield identical correlation energies. To the best of our knowledge, no such formal proof has been given before. 

As a method for calculating electronic excitation spectra, RPA requires the solution of 
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The matrices A, B, X, and Y are all ov x ov, where o and v are respectively the number of occupied and unoccupied 
spin-orbitals. The eigenvalue problem above can be completed by noting that if ( y* ) is an eigenvector with eigenvalue 

u>i, then ( ) is a l so an eigenvector, with eigenvalue — a>j. In the (real) canonical spin-orbital basis we use throughout 
this letter, we have 

Aiajb = (e a - £i)SijSab + ( ib || aj) , (2a) 
Biajb = (ij\\ab). (2b) 

Here, e p is a diagonal element of the Fock operator. Indices i, j, k, and I indicate occupied spin-orbitals, while a, 6, 
c, d indicate unoccupied spin-orbitals. For arbitrary spin-orbitals p, q, r, and s, the two-electron integral (pq || rs) is 
defined by 

(pq || rs) = (pq\ rs) - (pq\ sr) , (3a) 
(pq\rs) = J J dxi dx 2 p (xi) 9 (x 2 ) <?V( x i) <M*2) (3b) 

where x is a combined space and spin electron coordinate. 

The RPA correlation energy can be obtained by considering two harmonic excitation energy problems i 12 ' 15 RPA 
and the Tamm-Dancoff approximation (TDA) thereto, which sets B = and thus solves 

AZ = Zia (4) 

In the quantum chemistry community, TDA is also known as configuration interaction singles (CIS). While TDA 
includes only excitation operators, RPA also includes de-excitation operators which can be thought of as correlating 
the ground state. The ground state correlation energy in RPA is given by the difference between the zero point 
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energies of these two harmonic oscillator excitation problems with correlated (RPA) and uncorrelated (TDA) ground 
states. We thus have 

^r A = 5E'K-^) = ^Tr(a,-A). (5) 

i 

The prime on the summation means that we include only the positive excitation energies in defining uj. 

A different approach to calculating the correlation energy is coupled cluster theory. The simplest coupled cluster 
method includes only double excitations from the ground state, and is termed CCD. The CCD correlation energy is 

e c c cd = \Y,m ab ) ^ = lJ2( y i ab ) # . ( 6 ) 

where in the last equation we have used antisymmetry of t°*j under interchange of i with j or a with b. To determine 
the tfj, we solve the CCD equations in the spin-orbital basis (see, for example, Ref. [HI), 

= { ij \\ ab) + (e a + e b - e, - e^tfi + Jtg? ( ij \\ kl) + h% ( ab \\ cd) + \t$ { kl \\ cd) t% 

1 1 ( I \ ( ? ) 

- ^Pabtf ( kl || cd) ttf - -Biflft ( kl || cd) tH + l> ,!>.,/;, ({ ic || ah) + - ( kl || cd) #J , 

where internal indices (k, I, c, and d) are to be summed, and Py and P a b are permutation operators: (P a b9ac = Vac—Qbc, 
etc.). Keeping only particle-hole ring contractions, leads to what we shall here term "ring-CCD" (rCCD), 

= ( ij || ab) + tt (e c - e k ) 5 bc 5 jk + (e c - e k ) 5 ac 5 lk tf 3 + < ic \\ ak) t c k ) + t% < jc \\ bk) + t% ( kl |j cd) i*. (8) 

Defining tfj = Iio,j'6i an d using Eqn. [21 we obtain 

B + AT + TA + TBT = 0. (9) 

Removing the exchange integrals (i.e. setting (pq\\ rs) — > (pq\rs)) in Eqn. [5] gives us what we will call direct 
ring-CCD (drCCD), and in Eqn. [2] gives us direct RPA. Thus, Eqn. [9] holds both for rCCD and for direct rCCD with 
the A and B matrices defined as in RPA or direct RPA, respectively. In terms of B and T, the rCCD correlation 
energy is 

E r c CCD = ^Tr(BT), (10) 
while the drCCD correlation energy picks up an extra factor of two due to the different definition of B: 

E drCCD = I T r(BT). (11) 

We prove here that Eqn. [5] can be obtained from the RPA equations, and that with T thereby defined, the direct 
rCCD correlation energy of Eqn. [IT] is equal to the direct RPA correlation energy of Eqn. [5] 
We begin with the RPA equations, Eqn. Q] Multiplying on the right by X -1 , we have^i 



-B -a) (t) (t 



R, (12) 



where we have defined 

T = YX"\ (13a) 
R = XwX^. (13b) 

As seen below, T = YX _1 corresponds to the solution of Eqn. [9l Multiplying on the left by (T — 1) yields 

-b-a)(t) = < t (») 

Carrying out the matrix multiplications, we sec that this is just Eqn. [9l From Eqn. 1121 we have 

A + B T = R, (15) 
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whence 

Tr(B T) = Tr(R - A) = Tr(w - A). (16) 

The direct ring-CCD correlation energy is thus equal to the direct RPA correlation energy. The extra factor of 1/2 in 
the ring-CCD correlation energy on the right-hand-side of Eqn. [TO] makes the correlation energy exact to lead order, 
and it has been argued that it should therefore be included in defining the full RPA correlation energy. Sec Ref. N 
and references therein for discussion of this point. 

In order to obtain X and Y once we have T, we can use Eqn. [TO] to construct R. From Eqn. I13bl we can diagonalize 
R to get X. Once we have X and T, we simply use Y = T X to get Y. 

Direct RPA is commonly used in condensed matter physics, where the exchange terms are usually removed from the 
two-particle Hamiltonian (and treated as vertex corrections), and where typically scmilocal DFT orbitals and orbital 
energies (i.e. those coming from the local density approximation or a generalized gradient approximation) are used. 
The exchange-correlation energy in such a scheme is given by 

E xc = E« F +E d c RPA , (17) 

where E BF is the Hartree-Fock-type exchange energy with the semilocal orbitals and where "dRPA" indicates direct 
RPA. The pros and cons of keeping or neglecting vertex corrections in RPA correlation have been discussed in the 
literature £ 

Given that both ( ib | aj) and —tfj are positive definite for dRPA,— we can use Cholesky decomposition to write 

{ib\aj) = (ij\ab) = u? a uf b , (18a) 
-t$ = 0f a ef b , (18b) 
where A is to be summed. This leads to the drCCD equation (Eqn. [8] with no exchange integrals) becoming 



(n, A n, A o, A n, A B (i B A A , , B , , B _L (t A A n , B n , B d C d C \ 1 1 Q \ 

*ij = A~~afa V ia i h ~ ia kc kc J b ~ ia kc kc 3 b + ia kc kc ld ld 3 b > ^ ' 

ij 



with 



Defining 



Atf = + tj - e a - e b . (20) 

M AB =6lul, (21a) 
N AB =ule£ c , (21b) 
the construction of which scale as 0(ovc 2 ) where c = dim{A}, leads to 

^ = ^ K < b - ut n ab ef b e A m ab 9 B + e A u AB n bc e%) , (22) 

ij 

which can be solved by fixed point iteration with DIISi£ in 0(ovc 2 ) operations. Analytic energy gradients can also 
be carried out using the standard CC approach^ 

In the current framework, the cost of RPA is not much greater than that of MP2. The atomic orbital to molecular 
orbital integral transformation needed to build (ib\aj) scales as 0(N 5 ) for N atomic orbitals, and the Cholesky 
decomposition for dense (ib\aj) and tfj will scale worse than 0(N 4 ). However, transforming back into the atomic 
orbital basis (as in our AO-CC based formalismi^) will yield algorithms that scale near-linearly for sparse enough 
matrices^ 

The connection between the symplectic eigenvalue problem (Eqn. [1J and its associated Riccati equation (Eqn. [9]) 
is textbook material in Optimal Control Theory (see, for example, Ref. [2lh . Sanderson^ seems to have been the 
first to document this connection in the context of RPA; however, he neither mentions coupled cluster theory nor 
the agreement of correlation energies between RPA and rCCD. His assumption about commuting boson excitation 
operators leads to an RPA ground state representation that is correct only for two-electron systems^ 

In summary, we have offered an analytic proof that the excitation amplitudes of an approximate CCD model are 
related to the eigenvectors of the RPA model by T = YX" 1 , and that the ground state correlation energies of these 
two models are identical. This connection also lets us establish an 0(N i ) algorithm for the RPA correlation energy 
in a CC framework thanks to the mathematical properties of the solution (T < 0). 

This work was supported by the National Science Foundation (CHE-0807194 and CCF-0634902) and the Welch 
Foundation (C-0036). We thank Filipp Furche for providing benchmark numerical results of direct RPA correlation 
energies, and Gcorg Kresse for recently reviving our interest in this problem. 
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APPENDIX A: MATHEMATICAL DETAILS 

We here prove several statements about the solution of Eqn. [9l 



1. Symmetry of Y T X 



We begin by showing that X T Y = Y T X. Start with the RPA equation, Eqn. [TJ and multiply on the left by 
(Y T - X T ) to get 

Y T A X + X T A Y + Y T B Y + X T B X = (Y T X-X T Y)w. (Al) 
Since the left-hand-side is symmetric, we have 

(Y T X-X T Y)u; = w(X T Y- Y T X). (A2) 
Defining S = Y T X - X T Y, wc thus have 

Sw + t l >S = 0. (A3) 

In indicial form, this is 

S ij (u i +Uj) = 0. (A4) 
Since we have taken u>i positive, we must have S = 0, and hence X T Y = Y T X. 

2. Existence of X -1 

For positive definite B (true in direct RPA, but not in the full RPA), the existence of X -1 can be proven. Suppose 
that X z = for some vector z ^ 0. Multiplying both sides of the RPA equations by z would then give us 

BYz = Xwz, (A5a) 

-AYz = Yuz. (A5b) 

Since X T Y = Y T X, we would have 

z T Y T BYz = z T Y T Xwz=z T X T Y T wz = 0. (A6) 

Since B is positive definite, this implies that Yz = 0. But this would mean that (y) z = 0' contradicting the 
assumption that ( y ) is of full rank made in writing the eigenvalue problem. 

3. Symmetry of T 

Since X is nonsingular, and Y T X = X T Y, we have 

(X T )- 1 (Y T X-X T Y)X- 1 =0. (A7) 

Expanding the foregoing shows that 

T T T = 0. (A8) 

4. Negative Definiteness of T 

Since T is real and symmetric, we can diagonalize it with a unitary transformation U: T U = U A. Multiplying 
the drCCD equation on the left by a particular eigenvector and on the right by Ufe, we get 

dBU fc (l + A 2 fc )+2U'AU fc A fe = 0. (A9) 
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When B and A arc positive definite, as they are for direct RPA, we see that we must have Xk < for all k, and T is 
therefore negative definite. 
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